High-performance computing and modern numerical algorithms have made high-fidelity fluid-thermal analysis tractable in geometries of ever increasing complexity. Despite continued advances in these areas, direct numerical (DNS), large eddy simulation (LES), and even unsteady Reynolds-averaged Navier-Stokes (URANS) simulations of turbulent thermal transport remain too costly for routine analysis and design of thermal-hydraulic systems, where hundreds of cases must be considered. Reduced order models (ROMs) offer a promising alternative by leveraging expensive high-fidelity simulations (referred to as full order models or FOMs) to first extract a low-dimensional basis that captures the principal features of the underlying flow fields, and then construct computational models whose dimensions are orders of magnitude lower than the FOM dimension. In the numerical simulation of fluid flows, Galerkin ROMs (G-ROMs), which use data-driven basis functions in a Galerkin framework, have provided efficient and accurate approximations of laminar flows, such as the two-dimensional flow past a circular cylinder at low Reynolds numbers [1, 2]. However, turbulent flows are notoriously hard for the standard G-ROM. Indeed, to capture the complex dynamics, a large number [3] of ROM basis functions is required, which yields high-dimensional ROMs that cannot be used in realistic applications. Thus, computationally efficient, low-dimensional ROMs are used instead. Unfortunately, these ROMs are inaccurate since the ROM basis functions that were not used to build the G-ROM have an important role in dissipating the energy from the system [4].